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ABSTRACT 



u 

■ Reverberation mapping offers one of the best techniques for studying the inner 
|| regions of QSOs. It is based on cross-correlating continuum and emission-line light 

curves. New time-resolved optical surveys will produce well sampled light curves for 
many thousands of QSOs. We explore the potential of stacking samples to produce 
composite cross-correlations for groups of objects that have well sampled continuum 
. light curves, but only a few (~ 2) emission-line measurements. This technique exploits 

current and future wide-field optical monitoring surveys (e.g. Pan-STARRS, LSST) 
and the multiplexing capability of multi-object spectrographs (e.g. 2dF, Hectospec) to 
significantly reduce the observational expense of reverberation mapping, in particular 
at high redshift (0.5 to 2.5). 

We demonstrate the technique using simulated QSO light curves and explore the 
biases involved when stacking cross-correlations in some simplified situations. We show 
that stacked cross correlations have smaller amplitude peaks compared to well sampled 
correlation functions as the mean flux of the emission light curve is poorly constrained. 
However, the position of the peak remains intact. We find there can be 'kinks' in 
stacked correlation functions due to different measurements contributing to different 
parts of the correlation function. While the magnitude of the kinks must be fitted 
for, their positions and relative strengths are known from the spectroscopic sampling 
distribution of the QSOs making the bias a one-parameter effect. We also find the S/N 

■ in the correlation functions for the stacked and well-sampled cases are comparable for 
| the same number of continuum and emission line measurement pairs. 

& • Using the Pan-STARRS Medium-Deep Survey (MDS) as a template we show that 

cross-correlation lags should be measurable in a sample size of 500 QSOs that have 
weekly photometric monitoring and two spectroscopic observations. Finally we apply 
the technique to a small sample (42) of QSOs that have light curves from the MDS. 
We find no indication of a peak in the stacked cross-correlation. A larger spectroscopic 
sample is required to produce robust reverberation lags. 

Key words: quasars: emission lines, quasars: general, galaxies: nuclei, galaxies: ac- 
tive, galaxies: Seyfert 



1 INTRODUCTION 

The inner regions of active galactic nuclei (AGN) offer a 
unique opportunity to study matter within a few parsecs 
of a super-massive black hole. Reverberation mapping is 
designed to study (primarily) the broad-line region (BLR) 
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of AGN by measuring the inte raction between cont i nuum 
and broad-lin e flux variations (|Blandford fc McKeej Il982l : 
|Petersonlll993l '). The physical model assumes the BLR is pho- 
toionised by a UV continuum that is emitted from a much 
smaller radius. Variations in the ionising continuum produce 
equivalent variations in the broad emission line flux after a 
delay that can be associated with the light travel time. 

Reverberation mapping of a single system requires 
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many epochs of emission-line and continuum luminosity 
measures. A peak in the cross correlation between the two 
light curves indicates the time lag between continuum and 
emission-line variations. To date lags have be en measured for 
some tens of objec ts foll owing this approach JPeterson et al.l 
12004 iBentz et al.ll2006l . I2OI0I ; iDennev et aLll2006l . I2OI0F 

Reverberation mapping has led to significant advances 
in the u nderstanding of AGN (e.g. the radius- luminosity 
relation; IWandel et ail Il999l ; iKaspi et all 120001. stratifica- 
tion and kinematics of the BLR Peterson & Wandcl 1999, 
2000, black hole mass estimates IPeterson et alj |2004| etc.). 
However, these campaigns are observationally expensive and 
time consuming as they require many observations of indi- 
vidual objects over a long period of time. 

Rather than focusing on obtaining many hundreds of 
epochs of observations on single objects this paper explores 
an observational technique that is only now becoming possi- 
ble due to the new generation of time-resolved photometric 
surveys and the multiplexing capabilities of multiple-object 
spectrographs (MOS). 

Time resolved surveys, such as those bei ng performed 
with the Pan-STARRSl (PS1) telescope l|Kaiser et all 
2002), will measure the broad-band (continuum) light curves 
of many thousand of QSOs. The PS1 medium-deep survey 
(MDS) is taking images of ten 7deg 2 fields every few days 
in five photometric bands. Number counts im ply there are 
~ 50 QSOs with g < 22 in each of these fields (|Croom et al.l 
2009) making them easily surveyable with current MOSs 
in ~a night of observing time. Repeating the spectroscopic 
observations regularly would allow traditional reverberation 
mapping to be performed on samples of thousands of QSOs. 
However, this would require a large amount of observing 
time on highly subscribed telescopes to produce results. In 
this paper we look at the potential for reverberation map- 
ping in stacked samples of QSOs that only have a few (~ 2) 
spectroscopic epochs of data, but are coupled with well- 
sampled continuum light curves. 

In section [5] we outline the principle of stacking cross 
correlations, then in section [3] we simulate QSO light curves 
to illustrate the technique, in section [3] we present an em- 
pirical simulation of a QSO survey and the potential re- 
sults, in section [5] we apply the technique to early data 
from a spectroscopic survey of QSOs in the MDS region and 
in section [6] we summarise the results of our investigation. 
Throughout this paper we use a flat (fi m ,S7A) = (0.3,0.7), 
Ho ~ 70kms~ 1 Mpc -1 cosmology. 



2 OUTLINE OF APPROACH 

Observed light curves are made up of a series of (often un- 
evenly sampled) discrete measurements. Estimates of the 
cross-correlation between continuum and emission line light 
curves (denoted C and L below) are therefore limited by 
their sampling. A variety of techniques have been developed 
for estimating cross correlations with observational samples 
ilGaskell fc Sparkdl 19861 ; iGaskell fc Petersonll 19871 ; IZu et all 
l201ll ). In this pap er we will focus on the discrete cross co- 
variance function l|Edelson fc Krolik| [l988) as it lends itself 
simply to the stacking technique. The discrete cross covari- 
ance is calculated in terms of pairs of observations (d,Lj). 
Taking all pairs of continuum and emission line observations 



such that the time lag ti — tj is between r and r + St, then 
the cross covariance amplitude for lag r is estimated with 

*« — t j S[t,t + i5t] 

X{t)= Yl (Cl ~ g)(Lj ~ L) (1) 

. . Impair 

where nc and til are the number of continuum and emission 
line measurements respectively. St defines the bin size in 
the cross-corre l ation. Here we will use fixed bin sizes, but 
SCO lAlexanderl (|l997l ) for the use of variable bin sizes to 
optimise results. Equation[T]is the same as that for the cross 
correlation except that it is not normalised by the rms of 
each variable. St defines the bin size in the cross-correlation. 
In this paper we will focus on the covariance rather than the 
correlation function as they are almost equivalent in terms 
of the analysis presented here but using the cross covariance 
clarifies the discussion in later sections. 

The discrete cross covariance is a function on pairs of 
observations. Assuming all measurements have the same er- 
rors and well defined mean levels the variance of the cross 
covariance is inversely proportional to the number of data- 
data pairs. Therefore, desirable results require a large num- 
ber of pairs of observations at a wide range of time lags. 
Traditionally this is achieved through repeat observations 
of an object over a period of years to build up sufficient 
data at all time lags. However, there is no explicit reason 
that all of the data must come from a single source. Given a 
group of objects that have similar variability properties we 
are able to combine the data-data pairs from each object in 
the group to obtain a 'composite' cross covariance function 
with higher S/N than is obtainable for the individual ob- 
jects. Essentially this process is the same as calculating the 
covariance function for each object and then stacking them. 

This stacking technique has been used before for con- 
tinuum auto covariance/correlation analysis (lAlmaini et al.l 
120001 ; IVanden Berk et al.ll20o3 ; IWilhite et alj|2007l ). butto 
our knowledge it has not previously been applied to cross 
covariances and/or reverberation mapping. 

3 SIMULATING DATA 

To demonstrate the principle of stacking cross covariances 
we simulate QSO continuum and emission line light curves. 
We model the conti nuum as a first-order auto-regressive pro- 
cess following (e.g. iKellv et all 120091 ; iMacLeod et~al1 l20ld ; 
IZu et al.i r2012). Simulating this process requires a damping 
timescale r c and amplitude of short-time-scale variations a c 
that we fix for each simulation. 

The emission line flux at a given time is defined by 
convolving the continuum light curve w ith a given trans- 
fer fu nction that is a function of r (e.g. [Edelson fc Krolikl 
1988). The transfer function is not well constrained obser- 
vationally for AGN. For simplicity we use a Gaussian with 
unit area centred on r = t; with full-width-half-maximum 
n/2. Note that this is not a physically motivated transfer 
function, but is roughly similar to empirical transfer func- 
tions for Balmer lines that have been reverberation mapped 
l|Horne et al.lll99ll : IBentz et al.ll2010l ). 

After simulating the light curves we add random Gaus- 
sian noise into the continuum and emission line data with 
rms e c and e; to simulate uncorrelated effects and measure- 
ment error. 
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Table 1. Parameters used when simulating light curves (see text 
for the definitions of each parameter). In each case the parameter 
used in a simulation is drawn from a Gaussian distribution with 
mean and rms given in this table. Parameters that are required 
to be positive have their distribution functions truncated at zero. 



Parameter 




To 


n 


ec 


eL 


units 


flux 


days 


days 


flux 


flux 


mean 


5. 


20. 


40. 


0.2 


0.2 


rms 


1. 


5. 


5., 30. 


0.1 


0.1 



In total there are five input parameters. For each set 
of simulations we perform we will create many light curves. 
Each time we simulate a light curve we draw the input pa- 
rameters from Gaussian distributions with a mean and rms 
fixed for that set of simulations and table Q] gives a list of 
these parameters. 



3.1 Stacked vs. non-Stacked cross covariance 
functions 

To demonstrate the principle of stacking cross covariances 
we will concentrate on two simplified scenarios. In each we 
assume a QSO is observed photometrically (giving us contin- 
uum luminosity measures) every three days for six months 
of the year. In the first case (LI) we assume that at every 
date we get photometric data we also obtain a spectrum and 
emission-line flux measure. In the second case (L2) we only 
obtain line fluxes twice a year at either end of the contin- 
uum sampling. The LI scenario is designed to be equiva- 
lent to standard reverberation mapping with large numbers 
of emission line and continuum measurements, L2 demon- 
strates the stacking technique. We illustrate these situations 
in Fig.[T] The top panel in the figure shows the sampling over 
six months of the continuum and emission line in each of the 
cases. Below we show a simulated light curve of the contin- 
uum (solid line) and the emission line (dashed line) over the 
period (solid line) using the parameters from table [1] 

If we take just a single year's worth of data there are 61 
continuum measurements and 61 and 2 emission line mea- 
surements for LI and L2 respectively. In the case of LI we 
create 1000 simulated continuum and emission line light 
curves. We calculate the cross covariance function in each 
case and in Fig.[2Ja) we plot the average of the 1000 covari- 
ance functions along with the rms between simulations as 
error bars. 

Since there are ~ 30 times the number of spectroscopic 
measurements in LI we compare this with a case where 30 
QSOs have been sampled as in L2. The cross covariances 
from the 30 QSOs are then stacked to produce an ensemble 
covariance function. We simulate this situation 1000 times 
and plot the results in Fig. [2j6) and (c), the difference 
between these simulations is that the first has the rms of 
n — 5. days, while the second has t; = 30. days to illustrate 
the effect of increasing the scatter in reverberation lags. 

It is apparent that, in our simulations at least, we can 
reproduce the peak in the cross covariance function at t) = 
40 using stacked covariance functions. However, there are 
differences between the cross covariances for LI and L2. In 
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Figure 1. (bottom) Example of a simulated continuum (solid) 
and emission line (dashed) light curve. The dashed curve has been 
offset from the solid for clarity, (top) The dots in the top panel 
show how the continuum and emission line are sampled in the 
case of LI and L2. 



the following subsections we discuss the most prominent of 
these. 



3.1.1 Peak height in stacked vs. non-stacked cross 
covariances 

The peak in the covariance function for L2 is biased low be- 
cause we are under sampling the emission line light curve 
and so are unable to define the mean level precisely. The 
value of L in equation [T] is heavily biased by the two indi- 
vidual measurements of Lj, hence the covariance function 
is weaker. We can derive the expected level of this bias in 
the case of til measurements of the emission line luminosity. 
Taking the definition of the discrete cross covariance: 



T.L k 



hi k^j 



(d - C)(Lj 



= £■ 



-(Ci - C)(Lj - Lk^tj) 



(2) 



(note the sums here are over the same data points as equa- 
tion [1} L in equation [2] is now calculated over the ul — 1 
emission line measurements not including Lj and so is un- 
biased by Lj. Hence the amplitude of a cross covariance 
calculated from til emission line measurements will be pro- 
portional to (nij — 1) /ul ■ This is ~a factor of two in the case 
of LI compared to L2 in our simulations as is illustrated in 
Fig.pJ 

Note that had we calculated cross correlation rather 
than covariances this effect is not as important since the rms 
of the LiS (that normalises a cross correlation) is affected in 
the same manner as the covariance. 
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Figure 2. Comparison between the mean cross covariance functions for a single QSO with 61 epochs of spectroscopic data (LI; a) and 
30 QSOs with only two epochs (L2; b). (c) shows the same simulations as (b) except the distribution of lags used is 40±30days rather 
than ±5. In (d) the L2 simulation is re-run with rms(r;) = ±5. but the cross covariance functions are calculated using the known values 
of L and C rather than calculating them from the simulated observations. 
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Figure 3. The number of emission-line and continuum data- 
data pairs that are used in the cross-covariance calculation. We 
compare the LI and L2 examples (see text) with a bin size of six 
days. 



3.1.2 Errors in stacked vs. non-stacked cross covariances 

The rms error bars on the LI cross covariance are not all 
equal and decrease from left to right. In the case of L2 the 
error bars are roughly constant. The different ways that the 
emission line light curves are sampled for LI and L2 mean 
that the number of data-data pairs contribution to the co- 
variance function at a given t is not the same in each case. 

In Fig. we plot the number of data-data pairs con- 
tributing to a covariance function as a function of r for LI 
and L2. The constant time sampling of LI leads to a more 
concentrated distribution of lags. In the case of L2 the dis- 
tribution is flat while the total number of pairs is approxi- 
mately the same for LI and L2. This leads to an increase in 
the noise of LI cross covariances in the under-sampled areas 
of lag space around six months. 

In general the errors are considerably larger in the case 
of LI compared to L2. This is because LI over samples the 
light curve we have simulated and, while the number of mea- 
surements going into the LI and the stacked L2 cross covari- 
ances is roughly the same, the emission-line measurements 
are correlated in the case of LI, increasing the noise. 



3.1.3 The 'kink' in the L2 cross covariance 

The 'kink' in the L2 cross covariance occurs because we have 
only two measurements of the emission line luminosity, and 



because the mean of the continuum level is not accurately 
defined. 

At any lag in an individual cross covariance for L2 (be- 
fore stacking) there is only one data-data pair available to 
define the cross covariance. All of the positive lags use the 
first line measurement, and all of the negative lags use the 
second. The kink in Fig. b) occurs as we move from one 
regime to the other. To explain the kink we take the defini- 
tion of the discrete cross covariance in the case of just two 
line measurements: 

x{t) = y, (G-C) (£,--£) 

Impair 

= (C i -C){L j -(L 1 +L 2 )/2) (3) 

since there are only two emission line measurements and 
only one continuum — line pair contributes at a given lag. 
If we use Lji to denote the emission line measurement that 
isn't Lj then 

X{T) = Q.S{C i -C){L j -L jl ). (4) 

We now calculate the expectation value of the cross covari- 
ance. We use angled brackets to denote this expectation 
value that is averaged over a large number of realisations, 
this is distinct from the average of the individual contin- 
uum luminosities (C) that is calculated over the number of 
observations of a single object. 

(X(r)) = 0.5«C t Li) - (CiLj,) - (CL 3 ) + (CLj,)) 

we notate the covariance of two variables as cov(x,y) = 
(xy) - (x)(y) 

(X(t)) =0.5{{C i ){L j ) +cov(d,L i ) - {C l ){L ] ,)- cov(d, Lj,) 
-(C){Lj) -cov(C,Lj) + (C)(Lji) + cov(C, Lj,)) 

= 0.5(cov(Ci, Lj) — cov(Ci, Lj,) — cov(C, Lj) +cov(C, Lj,)). 

(5) 

The first term of equation [5] is the quantity we are try- 
ing to measure with the cross covariance. The other terms 
are biases. The factor of 0.5 is due to there only being two 
emission line measurements (see section I3.1.1[) . For the sec- 
ond term note that cov(C;, Li) < cov(C;, L2). However, this 
term will be small with respect to the other biases. For the 
other terms: cov(C, Li) will be small as C is calculated from 
the continuum light curve that is measured after L\ while L\ 
is dependent on the continuum level before it is measured. 
cov(C, 1/2 ) on the other hand is significant. This is partly 
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due to the manner in which we have constructed our simula- 
tion. The auto-regressive continuum light curve means that 
continuum points are covariant with their neighbours, as we 
may expect is the case in reality. Furthermore, the transfer 
function that is used to define the emission line luminosity 
makes L2 covariant with a number of the continuum points. 

Given the values in Table[T]we can calculate cov(C, Li), 
cov(C, £2), cov(Ci,Z/2) and cov(C nc ,L\) to get the ex- 
pected offset in the cross covariance between r < and 
r > 0. These are respectively 0.43, 5.04, 0.03, 0.0004 in 
the units used. Hence we expect the offset to be ~ 4.6 in 
Fig- [21 b). This is in good agreement with our simulations. 

In Fig.[2|d) we reproduce the simulations in (6), except 
that rather than calculating C and L from the simulated 
light curves, we use their correct values as defined in the 
simulation. In this case the bias is not apparent and there 
is no kink in the covariance function. 

3.2 Optimising the stacked results 

The major source of bias in the stacked covariance functions 
is the poorly defined mean levels, particularly L. The preci- 
sion of C can be improved simply with more measurements 
and, in the practical example we give below of the PAN- 
STARRS survey there will eventually be many more than 
the ~60 photometric measurements that we have assumed 
here. 

While more spectroscopic measurements would increase 
the precision of L this somewhat goes against the principle 
of the technique. Less biased values for L can be estimated 
from C and the global equivalent width distribution at fixed 
luminosity. However, this would also include a significant 
loss of precision and would smooth the kink in the correla- 
tion function at the expense of S/N. Since the time sampling 
of the continuum and emission line is known, the location 
of the kink and the shape (if not magnitude) of its effect 
on the covariance function can be estimated. Removing the 
bias from the results then requires a one-parameter fit to the 
data. In this case the magnitude of the bias can hold impor- 
tant information on the covariances between emission line 
and continuum emission in QSOs and, accurately measured, 
could help studies of transfer functions and the interactions 
between the accretion disk and BLR. 

The distribution of the reverberation lags of the quasars 
that are stacked in this manner has the effect of smoothing 
out the stacked covariance function. This broadens the peak 
in the stacked covariance function and hence reduces the 
signal in the peak. The two values of rras(ji) used in the 
above simulations are used to illustrate this effect with ex- 
treme values. The lack of large numbers of objects that have 
several reverberation mapped lines means that we cannot be 
certain of the distribution of lags for high redshift quasars. 
However, the H/3 line has been mapped for a significant num- 
ber of Seyferts. While there has been a considerable range of 
lags measured, almost all of this variation has been shown to 
be due to the radius-luminosity relation. iKaspi et al.l (120051 ) 
find only ~ 15 % intrinsic scatter around this relation. Fur- 
thermore, the small degree of scatter in Mgn and Civ line 
widths for the brightest quasars may indicate that there is 
even less intrin sic scatter in the radiu s-luminosity relation 
in that regime l|Fine et alj|2008t l2010t ). Hence it would be 
advisable to use quasars with a small range of luminosities 



when calculating a stacked cross covariance to improove the 
signal. Indeed this may be advisable anyway since one of the 
primary applications of this technique would be to evaluate 
the radius-luminosity relation for UV quasar lines. 

To get a better idea of what results we may expect from 
a feasible survey of quasars we design a simulation based on 
what is currently known about AGN. In the next section we 
present this simulation to give a realistic impression of the 
results that could be achieved with this technique. 



4 A PHYSICALLY MOTIVATED SIMULATION 
FOR STACKED REVERBERATION 
MAPPING WITH THE PAN-STARRS MDS 

In the remainder of this paper we make a preliminary ap- 
plication of the method in the Pan-STARRS medium deep 
survey (MDS). The MDS consists of ten 7deg 2 fields that 
are imaged every few nights in five photometric bands 
(gpirpiipizpiypi). The size of the fields, along with QSO 
number counts, make each field surveyable with the current 
generation of multi-object spectrographs in ~one night of 
observing time, yielding ~500 QSOs (with gp\ < 22). The 
regular photometric monitoring of the fields means that for 
every spectrum taken there are many continuum-emission 
line data pairs that can be used in cross covariance analy- 
ses. 

To evaluate the potential for reverberation mapping 
in the MDS fields we run simulations based on current 
knowledge of QSO variability parameters. We assume that a 
QSO survey is perfo rmed similar to the 2SLAQ QSO survey 
l|Croom et alj|2009f ) and simulate a single MDS field of data 
with the following prescription. 

• From 2SLAQ number counts we assume 500 QSOs to- 
tal. We therefore draw 500 QSOs at random from the 2SLAQ 
catalogue but only accept those objects at the right redshift 
to have Mgn in an optical spectrum (0.4 < z < 2.4). 

• For e a ch ob ject we use the BH ma ss estimates in 
iFine et all (|2008l ) and scaling relations from lMacLeod et al.l 
1 2010t) to obtain the continuum variability parameters ac 
and tc in the observed gpi-band assuming 0.3 and 0.15 dex 
scatter about the mean relations respectively. 

• We extrapolate from the i-band magnitude to 5100 A 
assuming a power-law continuum of fx oc A" 1 5 an d use 
the radius- luminosity relation from lKaspi et all |2005) with 
0.1 dex scatter to obtain the lag for the emission line re- 
sponse tl under t he assumption that the Mg 11 and H/3 lines 
have similar lags dMcLure fc Jarvis1l2002l ). 

• Timescales are converted to the observed frame and 
continuum light curves are simulated assuming 6 day sam- 
pling over an eight month period. Emission line flux values 
are calculated as previously with one spectrum at the be- 
ginning and end of the continuum points. 

• We then add random scatter to the continuum points 
based on the error on the 2SLAQ magnitude and add 10 % 
scatter to the emission line points to simulate measurement 
errors. 

In Fig. [4] we bin eight months worth of data by the 
ipi-band absolute magnitude and plot the stacked cross- 
covariances. The shaded areas show the regions within the 
rms of the mean cross-correlation when using 12 day bins 
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Figure 4. Simulated cross covariance functions between Mgn 
emission-line and gpi-band light curves for a survey of an MDS 
field. We assumed there would be 500 objects in the field and 
simulated individual QSO parameters by drawing randomly from 
the 2SLAQ catalogue. The objects are split into four bins by 
absolute magnitude to show the tendency towards longer lags in 
brighter objects. 
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Figure 5. Simulated cross-covariance function for a single mag- 
nitude bin in one of our simulations (points) and the functional 
fit to it (solid line). Here a Gaussian fits the simulated points and 
the kink at r = days is corrected by the step function. 



in the discrete cross-covariance. The figure shows that the 
tendency towards longer lags in brighter QSOs may be de- 
tectable in a single year's worth of data on a single MDS 
field given the above assumptions. 

For each individual realisation in our simulations we try 
to find the peak in the cross-covariance functions. We fit a 
Gaussian, plus an offset, plus a step at r = Odays. The step 
function corrects for the kink in the covariance functions at 
t — Odays. In Fig. [5] we show the simulated cross-covariance 
for a single magnitude bin in one of our realisations along 
with the fitted Gaussian+step. 

Using the fit Gaussians we find the location of the peak 
in the cross correlations. The distributions of these values 
are shown in Fig. [5] for each of the magnitude bins. Note 
that we do not constrain any of the parameters in the fits 
but still find a peak in the range of timescales sampled for 



Figure 6. The distribution of the peaks in individual simulated 
cross-covariances, binned by ipi-band magnitude, as defined by 
Gaussian fits, the tendency for longer lags in the brighter bins 
is clear here. The brightest bin has lags that are too long to be 
properly sampled and we only find a good peak in ~ 40 % of the 
simulations. 



almost all (>99 %) of the realisations except in the bright- 
est magnitude bin. In the brightest bin the timescales are 
longer and we only obtain a lag in ~ 40 % of the realisations. 
Note that in some cases a Gaussian gives a poor automated 
fit to the cross-covariance. More careful fitting in individ- 
ual cases would give better results. The mean tl values in 
the simulated light curves are 27, 54, 97 and 227 days for 
the faintest to brightest magnitude bins. The fact that we 
measure considerably smaller lags than this is due to the 
fainter objects in each bin being more variable and hence 
contributing more to the stacked correlation function. 



4.1 The potential of stacked cross-covariance 
functions 

We have demonstrated that, at least in our simulations, it 
is possible to retrieve an average radius-luminosity relation 
for a sample of quasars. If possible this would be of con- 
siderable interest for quasar studies. The radius-luminosity 
relation has not been strongly constrained for any emis- 
sion lines other than H/3. However, UV lines (notably Mgn 
and C iv) are commonly used to es timate virial black hole 
masses in high-red shift quasars (e.g. lMcLure fc Jarvi 3 120021 ; 
IVester gaard 2002]). These virial mass estimates are based on 
the assumption that a radius-luminosity relationship exists 
for these UV lines equivalent to that measured for H/3. An 
observational determination of these relations, even if only 
in averaged stacks, would be of great use in determining 
black hole masses in high-redshift quasars. 

The simulations performed in section U assumed a par- 
ent sample of 500 quasars. The MDS survey as a whole would 
contain ~ 5000 quasars to the flux limit we assume. It is not 
beyond current telescopes to survey this number of objects 
~ twice yearly. Such a dataset, in particular if built up over 
several years, would offer the potential to make extremely 
high- precision stacked covariance functions. These could of- 
fer a unique opportunity to study the transfer function and 
hence the structure of the BLR although such studies would 
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have to be mindful that convolved into the covariance func- 
tion is the distribution of different lags and transfer func- 
tions of the different quasars in a stack. 

Multiple lines could be mapped for the same stacks of 
quasars. These data give relative information on the stratifi- 
cation of the BLR and where the various emission zones are 
in quasars. Combined with dynamical measurements (i.e. 
line profiles) these give indications of the dominant motions 
in the BLR. 

The potential gains from being able to reverberation- 
map high-z quasars are significant. Above we have outlined 
just a few. This paper presents a technique that we believe 
offers a feasible rout towards obtaining these results. 

Recent studies have also suggested either narrow- or 
broad-band 'photometric' reverberation mapping as an 
observationally cheaper m ethod for obtaining reverbera- 
tion lags for quasars (e.g. ( Cherepashchuk fc Lvutvil 1 19731 ; 
lHaas et aljkulll ; IChelouche fc Danidl2012ft ). While neither 
technique has been proven for large-scale samples the rela- 
tive gains of each technique are somewhat unclear. Narrow- 
band photometric mapping requires the quasars are at the 
correct redshift for the filter, and so cannot be applied on 
such large scales as the other two methods. Broad-band pho- 
tometric reverberation mapping is more efficient than the 
technique described here in terms of the spectroscopy re- 
quired. However, a single epoch of spectroscopy is required 
to obtain accurate redshifts. The complexity of decoupling 
continuum and emission line light curves and the effects of 
having several lines in a single filter requires extremely accu- 
rate photometry and complex reduction techniques. On the 
other hand broad-band photometric reverberation mapping 
does not require the complexity of obtaining flux-calibrated 
fibre spectra with high enough S/N to calculate precise line 
fluxes. Finally, a major gain for stacked covariance functions 
is the flat sampling distribution (Fig. [3} • Bright quasars can 
have continuum timescales of years, and correspondingly 
large BLR lags. Obtaining good constraints on these lags re- 
quires sampling many times the continuum+lag timescale in 
classical reverberation mapping. However, in our technique 
one needs only sample 1 x the continuum+lag timescale and 
then information is built up through stacking. 

The problem of obtaining robust reverberation mapping 
results in the high-redshift Universe is a problem that may 
be solved with current and planned time-resolved photomet- 
ric surveys. In the rest of the paper we apply our technique 
to a small number of objects that already have the correct 
observations available. 



5 PAN-STARRS MDS EARLY DATA 

In this section we derive MDS light curves for a sample of 
QSOs that have > 1 spectra, and calculate the stacked cross- 
covariance function for the sample. 




rms (Mag.) 

Figure 7. The rms of the calibration offset between PS1 gpi and 
SDSS g-band photometry in the MD03 field. In cases where the 
rms is > 0.1 Mag the calibration is deemed to be suspect and the 
skyccll is rejected. This corresponds to ~ 5 % of the individual 
skycells. 



photometric bands while not affected by the sun or moon. 
On each night eight dithered exposures are taken and com- 
bined to form a n ightly stack (see e.g. iKaiser et alj 120101 ; 
iTonrv et al.l 120121 for further details of the PS1 telescope, 
observing strategy and data processing). We used psphot, 
part of the stan dard PS1 Image Processing Pipeline system 
(Magnicr 2006), to extract point-spread- function photom- 
etry from nightly stacked images of the MDS fields. Each 
nightly stack is divided into ~ 70 skycells. W e calibrate each 
of th e se separately usin g SDSS photometry (|Fukugita et al.l 
1 19961 ; lYork et aljfeoOCh of moderately bright (16 < Mag. < 
18.5) point sources (defined as having type = 6 in the 
SDSS database). The bright magnitude cut is used since 
the brightest objects can create artifacts in the PS1 images. 
We use 18.5 as the faint cut so that we have a large number 
of objects used in the calibration of each skycell while en- 
suring high-precision photometry for each object (95 % have 
errors< 0.01 Mag.). 

We found that in some individual cases the flux calibra- 
tion was unstable due primarily to artifacts around objects 
used in the flux calibration in the PS1 imaging. To calibrate 
our data we measure an average offset between the PS1 and 
SDSS photometry, we also record the rms around this aver- 
age for each skycell we calibrate. Fig.Qshows a histogram of 
the rms of the calibration offset for each epoch of each sky- 
cell in the gpi-band for the MD03 field. The tail to larger 
rms' is indicative of PS1 skycells that have poor magnitudes 
for the stars used in the flux calibration. We only accept 
skycells with rms < 0.1 Mag. This cut removes ~ 5 % of the 
stacked skycells from our sample. 



5.1 PS1 light curves 

The PS1 telescope (jHodapp et al.|[2004 ) is performing a se- 
ries of time-resolved photometric surveys of the northern 
sky. We are particularly interested in the MDS as offering 
the best opportunity for our analysis. Images of the ten 
MDS fields are taken every four-five nights in each of the 



5.2 Hectospectra 

Spectra of QSOs in the MDS fields are being taken as part 
of an ongoing project to study the variability of QSOs. QSO 
candidates are selected for spectroscopy using current pho- 
tometric databases of QSOs based o n SDSS photometry 
l|Richards et al.ll2009l ; lBovv et al.ll201lft in addition to which 
point sources that correspond to X-ray sources, variability 
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selected objects and UVX selected objects are targeted for 
spectroscopy. 

The MDS fields are being surveyed with the Hectospec 
instrument on the MMT. Each MDS field is tiled with seven 
MMT pointings. Exposures are ~1.5h in length meaning 
that an MDS field (~ 500 QSOs) can be surveyed in ~ 
f night of on-sky observing time. 

The spectra are extracted and red uced using standard 
Hectospec pipelines (|Mink et alj|2007h . They are then flux 
calibrated using observations of F stars in the same fields. 
These stars are compared with a grid of model stellar spec- 
tra, created using the spectral synthes is code spectrum 



jGrav fc Corballvl 19941; iGrav et al.ll200ll ) with models from 
ICastelli fc Kuruc 5 |2004i r to correct for the response of 
the Hectospec instrument. Absolute flux calibration is then 
made by fitting to the r-band SDSS magnitudes of the stars. 
Errors in the response correction are typically < 10 % over 
the main part of the spectrum (~4000-8500 A). However, 
uncertainty on the absolute flux calibration can have a larger 
effect on the line fluxes we measure. This depends strongly 
on the number of stars used in the calibration and is typi- 
cally < 10 %. However, in a handful of Hectospec fields some 
stars are significantly off the average calibration. In general 
the spectra of these stars are fainter than expected and may 
be affected by small positioning errors. These stars are re- 
moved manually from the calibration. However, the same 
problem may affect QSO spectra in our sample. 

5.3 Our sample 

In total spectra of 855 (g < 22) QSOs have been obtained, 
primarily in MDS fields MD03 and MD07. More informa- 
tion on the selection of this QSO sample and the spectra 
will be given in an upcoming paper. Further QSO spec- 
tra were kindly taken by the MDS transient team. Most 
QSOs in the sample have only one spectrum and so cannot 
be used to cross correlate with the continuum observations. 
However, 82 of these objects ar e in the SDSS DR7 QSO 
catalogue ({Schneider et al.ll2010l ). These objects have spec- 
tra taken with the Sloan telescope giving us two spectra 
over a time baseline of ~ 10 years (for d etails of the Sloan 
spectrograph and SPS S QSO selection see | Gunn et al.l 2006: 
IStoughton et all |2002| ; iRichards et all 120021 ). Furthermore, 
we have taken > 1 spectrum of 59 QSOs in our sample giv- 
ing us a spectroscopic baseline from ~50 to ~100days. In 
Fig. [8] we show the rpi-band light curve for 4 QSOs selected 
from our sample as examples. 
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Figure 8. Four example rpi-band light curves of QSOs in our 
sample. The top two and bottom two come from different MDS 
fields (MD03 and MD07 respectively). All of the objects show 
evidence for variability over the ~ 1 year of monitoring. Dashed 
lines in each plot indicate the SDSS r-band magnitude. 




5.4 Emission line fluxes 

In Fig. [5] we show the redshift distribution of the 138 QSOs 
with > 1 spectrum. Their distribution is relatively typical of 
that of the SDSS (and other optically selected) QSO cata- 
logue with the vast majority between z ~ 0.5 and 2. This 
corresponds roughly to the range of redshifts over which the 
Mgn line is redshifted into optical spectra and we focus on 
the Mg n line for the rest of this work. 

We fit th e Mg H line following the prescription outlined 
in iFine et all |2008l ) and each fit is manually inspected to 
check the reliability. We calculate emission line flux and er- 
ror dir ectly from th e cont inuum-subtracted spectrum fol- 
lowing [Cardiereinal| ||l998h . The typical S/N of these lines 



z 

Figure 9. The redshift distribution of the 138 objects that have 
MDS light curves and > 1 spectroscopic observation from the 
MMT and SDSS. The distribution is typical of optically/UV se- 
lected samples and is highly incomplete in the interval 2 < z < 3. 

is ~ 3-30 for the SDSS spectra and ~ 10 - 200 for 
the Hectospec spectra. Hence in the case of the Hectospec 
spectra the error on the line fluxes are dominated by the 
flux-calibration errors rather than the spectral S/N. The 
SDSS spectra, that have a flux calibration error of ~ 5 % 
l|Adelman-McCarthv et al"]|2008r ). are more typically domi- 
nated by the statistical noise in the spectra. 
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Figure 10. This figure compares the Mgn line flux measured 
from the archive SDSS spectra and during the Hectospec survey. 
Errors on the SDSS measurements are dominated by the spectral 
S/N, while in the Hectospec observations they are dominated by 
uncertainty in the flux calibration. 

Through visual inspection we remove a number of spec- 
tra that have i) broad-absorption lines, ii) unusual spectra 
that are poorly characterised by our fitting or iii) are affected 
by residual sky/telluric features, fn Fig.[lO]we show the mea- 
sured SDSS and Hectospec fluxes for the 42 QSOs that had 
good Mgn flux measurements from each spectrum. These 
span ~two orders of magnitude in luminosity and cover the 
redshift range 0.45 < z < 1.68. 

5.5 The binned cross covariance 

We calculate the discrete covariance function for the 42 ob- 
jects with two good Mg n flux measurements. In the calcula- 
tion we convert the emission-line fluxes to magnitudes, hence 
we are measuring the fractional (rather than absolute) cross 
covariance. We perform the calculation in bins with width 
20 days in the rest frame of the observed QSOs. In Fig. [TJJ 
we show the covariance function along with the number of 
data-data pairs in each time-lag bin. The small number of 
objects and spectra means that there is little signal in the 
cross covariance. Furthermore, due to the observing times 
we have a gap in our time sampling at ~ 90 days. 

The small number of objects involved means that we 
were unlikely to find a lag in this data. Fig. [TJJ is given 
rather as an example of what can be done with current data. 
Note that we do not expect the type of step bias shown in 
Fig. [2jc) since, in those objects that have spectra from the 
SDSS, only one of the two spectra is contributing to the 
correlation function at these lags. 



6 SUMMARY 

We have demonstrated a technique for stacking cross- 
correlations and covariances. We focus on QSO emission line 
and continuum light curves in the case when only one is well 
sampled. We demonstrate the technique via a suite of simple 
simulations that highlight some of its biases and limitations 
as well as its advantages. While the stacked analyses show 
smaller peaks and can produce erroneous steps in the results 




200 



t = t c — t L (rest frame days) 



Figure 11. The stacked discrete cross covariance function for the 
42 objects in our final sample. The bottom panel shows the num- 
ber of emission-line and continuum data-data pairs contributing 
to each bin. The bins are 20 days (rest-frame) in width. 

due to the limited time-sampling, the position of the peak 
remains unbiased. Furthermore, given similar total numbers 
of emission line and continuum measurements the technique 
gives comparable S/N compared to the classic, unstacked 
approach. 

We focus on the PAN-STARRS MDS and show that 
it may be possible to measure average reverberation lags 
for QSOs in these fields with a relatively small investment 
of telescope time based on empirical simulations of QSOs in 
these fields. Finally we performed a stacked cross-covariance 
analysis on 42 QSOs from the MDS that have well sampled 
continuum emission and > 2 spectra from our observations 
and the SDSS archive. We find no indication of a peak in 
these data although the small numbers mean this would have 
been unlikely. In the near future multi-epoch spectroscopic 
observations of MDS fields would be required to allow for 
stacked lags to be measured. 

We note that the same technique could be applied 
outside of reverberation mapping. The relationship be- 
tween X-ray and opt i cal variability is complex in QSOs 
JShemmer et alj 120031 ; iMarshall et al. I l2008l ; lArevalo et al.l 
2009) and lags between X-ray and optical light curves can be 
used to investigate the dominant emission processes at these 
wavelengths. Here the highly sampled optical data could be 
combined with wide-field X-ray imaging data (e.g. XMM) as 
an efficient means of producing optical to X-ray cross corre- 
lations. The gain in efficiency would be proportional to the 
number of X-ray QSOs that can be simultaneously imaged. 
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